/* "Efficiency and water use: Dynamic effects of irrigation technology adoption"
by Micah Cameron-Harp and Nathan Hendricks

Code written by Micah Cameron-Harp
May 16th, 2024

This do file reads the TWFE and Callaway & Sant'Anna dynamic treatment
	effect estimates and generates figures 4 and 5
*/

/* NOTE - Global macros for directories are created in the main_text_results.do 
	file which calls this individual .do file */
	
*Create log file
log using "$dr_output_log\figures4and5", replace

/***************************** Flood to CP/LEPA *******************************/
/* Import the flood to cp or lepa event study coefficients */
	*Start with TWFE
	import excel using "$dr_temp\twfe_event_study.xlsx", firstrow sheet(floodtocplepa_91_15_nyt) clear
	gen estimator = "twfe"
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	rename (ci_low_95 ci_high_95) (lb_estimate ub_estimate)
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep dep_var estimator effect_estimated scaled_estimate scaled_lb scaled_ub estimate lb_estimate ub_estimate
	save "${dr_temp}\floodtocplepa_es_twfe.dta", replace

	*Now do CS
	import excel "${dr_temp}/cs_floodtocporlepa_91_15_nyt.xlsx", firstrow sheet(cs_event) clear
	gen estimator = "c&s"
	foreach var of varlist effect_estimated estimate se_estimate c_crit_val panel_mean_dep_var {
		destring `var', replace
	}
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	rename (ci_95_lb ci_95_ub) (lb_estimate ub_estimate)
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep dep_var estimator effect_estimated scaled_estimate scaled_lb scaled_ub estimate lb_estimate ub_estimate
	save "${dr_temp}\floodtocplepa_es_cs.dta", replace

	*Append the TWFE results 
	append using "${dr_temp}\floodtocplepa_es_twfe.dta"
	
	*Generate column indicating tech transition and save
	gen transition = "flood_to_cp"
	save "${dr_temp}\floodtocplepa_es_twfecs.dta", replace 
 
********************************************************************************
/* Figure 4 - Plot dynamic treatment effect estimates for flood to center-pivot
or LEPA conversion in percent terms for TWFE and CS estimators */
********************************************************************************
*Sort 
sort dep_var effect_estimated

*Now offset the location of scaled_estimates for C&S so they don't overlap
replace effect_estimated = effect_estimated + 0.5 if estimator=="c&s"

*Replace extreme ub and lb estimates with -50 and +50 to make graphs readable 
gen extreme = scaled_ub > 40 | scaled_lb < -40
replace scaled_ub = 40 if scaled_ub > 40
replace scaled_lb = -40 if scaled_lb < -40

/* Graph C&S, C&D, and TWFE together */
*AF_USED_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="af_used" & extreme!=1, msize(.3) lwidth(.1) lcolor(orange%75)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="af_used", lwidth(.1) lcolor(orange%75) msize(.4) mcolor(orange%75)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="af_used",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%75)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="af_used" & extreme!=1,  msize(.3) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="af_used", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="af_used",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Flood to center pivot", size(small)) ///
			subtitle("Acre-feet withdrawn", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-24 -20(5)20 23, nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(5 2) label(5 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
			graphregion(color(white) margin(zero)) name(es_af_used, replace)

*ACRES_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="acres_irr" & extreme!=1, msize(.3) lwidth(.1) lcolor(orange%75)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="acres_irr", lwidth(.1) lcolor(orange%75) msize(.4) mcolor(orange%75)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="acres_irr",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%75)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="acres_irr" & extreme!=1,  msize(.3) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="acres_irr", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100)  msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="acres_irr",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Flood to center pivot", size(small)) ///
			subtitle("Acres irrigated", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-24 -20(5)20 23, nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(5 2) label(5 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							    label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_acres_irr, replace)

*DEPTH APPLIED
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="depth_applied" & extreme!=1, msize(.3) lwidth(.1) lcolor(orange%75)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="depth_applied", lwidth(.1) lcolor(orange%75) msize(.4) mcolor(orange%75) msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="depth_applied",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%75)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="depth_applied" & extreme!=1,  msize(.3) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="depth_applied", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="depth_applied",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Flood to center pivot", size(small)) ///
			subtitle("Depth applied", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-24 -20(5)20 23, nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(5 2) label(5 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_depth_applied, replace)
			
/* Combine the graphs */
grc1leg2 es_af_used es_acres_irr es_depth_applied, ///
	rows(3) cols(1) legendfrom(es_af_used) ///
	maintotoptitle maintitlefrom(es_af_used) ///
	xtob1title xtitlefrom(es_af_used) ///
	ytol1title ytitlefrom(es_af_used) ///
	labsize(vsmall) graphregion(color(white)) xsize(6.5) ysize(8) iscale(.8) ring(1)
graph export "$dr_output_main\figure4.tif", width(9750) height(11850) replace

/***************************** CP to LEPA *******************************/			
/* Import the cp to lepa event study coefficients */
	*Start with TWFE
	import excel using "$dr_temp\twfe_event_study.xlsx", firstrow sheet(cptolepa_91_19) clear
	gen estimator = "twfe"
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	rename (ci_low_95 ci_high_95) (lb_estimate ub_estimate)
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep dep_var estimator effect_estimated scaled_estimate scaled_lb scaled_ub estimate lb_estimate ub_estimate
	save "${dr_temp}\cptolepa_es_twfe.dta", replace

	*Now do CS
	import excel "${dr_temp}/cs_cptolepa.xlsx", firstrow sheet(cs_event) clear
	gen estimator = "c&s"
	foreach var of varlist effect_estimated estimate se_estimate c_crit_val panel_mean_dep_var {
		destring `var', replace
	}
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	rename (ci_95_lb ci_95_ub) (lb_estimate ub_estimate)
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep dep_var estimator effect_estimated scaled_estimate scaled_lb scaled_ub estimate lb_estimate ub_estimate
	save "${dr_temp}\cptolepa_es_cs.dta", replace

	*Append the TWFE  results 
	append using "${dr_temp}\cptolepa_es_twfe.dta"	
		
	*Generate column indicating tech transition and save
	gen transition = "cp_to_lepa"
	save "${dr_temp}\cptolepa_es_twfecs.dta", replace 

********************************************************************************
/* Figure 5 - Plot dynamic treatment effect estimates for center-pivot to LEPA
conversion in percent terms for TWFE and CS estimators */
********************************************************************************
*Sort 
sort dep_var effect_estimated

*Now offset the location of scaled_estimates for C&S so they don't overlap
replace effect_estimated = effect_estimated + 0.5 if estimator=="c&s"

*Replace extreme ub and lb estimates with -50 and +50 to make graphs readable 
gen extreme = scaled_ub > 40 | scaled_lb < -40
replace scaled_ub = 40 if scaled_ub > 40
replace scaled_lb = -40 if scaled_lb < -40

/* Graph C&S together with C&D */
*AF_USED_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="af_used" & extreme!=1, msize(.3) lwidth(.1) lcolor(orange%75)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="af_used", lwidth(.1) lcolor(orange%75) msize(.4) mcolor(orange%75)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="af_used",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%75)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="af_used" & extreme!=1,  msize(.3) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="af_used", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="af_used",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Traditional center pivot to LEPA", size(small)) ///
			subtitle("Acre-feet withdrawn", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-28 -20(5)20 27,  nogrid) ///
			xscale(range(-28.5 27.5)) ///
			legend(order(8 5 2) label(5 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
			graphregion(color(white) margin(zero)) name(es_af_used, replace)

*ACRES_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="acres_irr" & extreme!=1, msize(.3) lwidth(.1) lcolor(orange%75)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="acres_irr", lwidth(.1) lcolor(orange%75) msize(.4) mcolor(orange%75)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="acres_irr",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%75)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="acres_irr" & extreme!=1,  msize(.3) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="acres_irr", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="acres_irr",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Traditional center pivot to LEPA", size(small)) ///
			subtitle("Acres irrigated", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-28 -20(5)20 27,  nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(8 5 2) label(5 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_acres_irr, replace)

*DEPTH APPLIED
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="depth_applied" & extreme!=1, msize(.3) lwidth(.1) lcolor(orange%75)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="depth_applied", lwidth(.1) lcolor(orange%75) msize(.4) mcolor(orange%75)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="depth_applied",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%75)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="depth_applied" & extreme!=1,  msize(.3) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="depth_applied", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="depth_applied",  msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Traditional center pivot to LEPA", size(small)) ///
			subtitle("Depth applied", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-28 -20(5)20 27,  nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(8 5 2) label(5 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_depth_applied, replace)
			
/* Combine the graphs */
grc1leg2 es_af_used es_acres_irr es_depth_applied, ///
	rows(3) cols(1) legendfrom(es_af_used) ///
	maintotoptitle maintitlefrom(es_af_used) ///
	xtob1title xtitlefrom(es_af_used) ///
	ytol1title ytitlefrom(es_af_used) ///
	labsize(vsmall) graphregion(color(white)) xsize(6.5) ysize(8) iscale(.8) ring(1)
graph export "${dr_output_main}\figure5.tif", width(9750) height(11850) replace

*Close log
log close
